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Abstract 

To obtain the initial pressnre from the collected data on a planar sensor 
arrangement in Photoacoustic tomography, there exists an exact analytic 
frequency domain reconstruction formula. An efficient realization of this 
formula needs to cope with the evaluation of the dataaAZs Fourier trans¬ 
form on a non-equispaced mesh. In this paper, we use the non-uniform 
fast Fourier transform to handle this issue and show its feasibility in 3D 
experiments. This is done in comparison to the standard approach that 
uses polynomial interpolation. Moreover, we investigate the effect and 
the utility of flexible sensor location on the quality of photoacoustic im¬ 
age reconstruction. The computational realization is accomplished by 
the use of a multi-dimensional non-uniform fast Fourier algorithm, where 
non-uniform data sampling is performed both in frequency and spatial do¬ 
main. We show that with appropriate sampling the imaging quality can 
be significantly improved. Reconstructions with synthetic and real data 
show the superiority of this method. Keywords: Image reconstruction, 
Photoacoustics, non-uniform FFT 


1 Introduction 

Photoacoustic tomography is an emerging imaging technique that combines the 
good contrast of optical absorption with the resolution of ultrasound images 
(see for instance dH). In experiments an object is irradiated by a short-pulsed 
laser beam. Depending on the absorption properties of the material, some light 


energy is absorbed and converted into heat. This leads to a thermoelastic ex¬ 
pansion, which causes a pressure rise, resulting in an ultrasonic wave, called 
photoacoustic signal. The signal is detected by an array of ultrasound trans¬ 
ducers outside the object. Using this signal the pressure distribution at the time 
of the laser excitation is reconstructed, offering a 3D image proportional to the 
amount of absorbed energy at each position. This is the imaging parameter of 
Photoacoustics. Common measurement setups rely on small ultrasound sensors, 
which are arranged uniformly along simple geometries, such as planes, spheres, 
or cylinders covering the specimen of interest (see for instance [2I1E01I221II911I] 
below). 

For the planar arrangement of point-like detectors there exist several ap¬ 
proaches for reconstruction, including numerical algorithms based on filtered 
backprojection formulas and time-reversal algorithms (see for instance [20l fTTl 

The suggested algorithm in the present work realizes a Fourier inversion for¬ 
mula (see 0 below) using the non-uniform fast Fourier transform (NUFFT). 
This method has been designed for evaluation of Fourier transforms at non- 
equispaced points in frequency domain, or non-equispaced data points in spa¬ 
tial, respectively temporal domain. The prior is called NER-NUFFT (non- 
equispaced range-non uniform FFT), whereas the latter is called NED-NUFFT 
(non-equispaced data-non uniform FFT). Both algorithms have been introduced 
in [7] . Both NUFFT methods haven proven to achieve high accuracy and simul¬ 
taneously reach the computational efficiency of conventional FFT computations 
on regular grids [7]. This work investigates photoacoustic reconstructions from 
ultrasound signals recorded at non-equispaced positions on a planar surface. To 
the best of our knowledge, this is a novel research question in Photoacoustics, 
where the use of regular grids is the common choice. For the reconstruction 
we propose a novel combination of NED- and NER-NUFFT, which we call 
NEDNER-NUFFT, based on the following considerations: 

1. The discretization of the analytic inversion formulas (see 0 ) contains 
evaluation at non-equidistant sample points in frequency domain. 

2. In addition, and this comes from the motivation of this paper, we consider 
evaluation at non-uniform sampling points. 

The first issue can be solved by a NER-NUFFT implementation: For 2D pho¬ 
toacoustic inversion with uniformly placed sensors on a measurement line such 
an implementation has been considered in . This method was used for biolog¬ 
ical photoacoustic imaging in m- In both papers the imaging could be realized 
in 2D because integrating line detectors mm have been used for data record¬ 
ing. In this paper, however, the focus is on 3D imaging, because measurements 
are taken with point sensors. Experimentally, we show the applicability and 
superiority of the (NED)NER-NUFFT reconstruction formula in three spatial 
dimensions, compared with standard interpolation FFT reconstruction. To be 
precise, in this paper we conduct three dimensional imaging implemented us¬ 
ing NEDNER-NUFFT, with ultrasound detectors aligned non-uniformly on a 


measurement plane. To easily assess the effect of a given arrangement, also 2D 
numerical simulations have been conducted, to support the argumentation. We 
quantitatively compare the results with other computational imaging methods: 
As a reference we use the fc-wave toolbox HZI with a standard FFT implementa¬ 
tion of the inversion algorithm. The NEDNER-NUFFT yields an improvement 
of the lateral and axial resolution (the latter even by a factor two). 

The outline of this work is as follows: 

In Section we outline the basics of the Eourier reconstruction approach by 
presenting the underlying Photoacoustic model. We state the Eourier domain 
reconstruction formula 0 in a continuous setting. Moreover, we figure out two 
options for its discretization. We point out the necessity of a fast and accurate 
algorithm for computing the occurring discrete Eourier transforms with non- 
uniform sampling points. 

In Section we briefly explain the idea behind the NUFFT. We state the 
NER-NUFFT (subsection 3.1) and NED-NUFFT (subsection 3.2) formulas in 
the form we need it to realize the reconstruction on a non-equispaced grid. 

In Section we discuss the 3D experimental setup. The NER-NUFFT is 
compared with conventional FFT reconstruction. A test chart is used to quantify 
resolution improvements in comparison to the k-wave FFT reconstruction with 
linear interpolation. In axial direction this improvement was about 170% while 
reducing the reconstruction time by roughly 35%. 

In Section [5] we then turn to the NEDNER-NUFFT in 2D with simulated 
data, in order to test different sensor arrangements in an easily controllable 
environment. An equiangular arrangement turns out to yield an improvement 
of over 40% compared to the best choice of equispaced sensor arrangement. 
Furthermore, we use the insights gained from the 2D simulations to develop 
an equi-steradian sensor arrangement for our 3D measurements. We apply our 
NEDNER-NUFFT approach on these data and quantitatively compare the out¬ 
comes with reconstructions from equispaced data obtained by the NER-NUFFT 
approach. Our results show a significant improvement of the already superior 
NER-NUFFT. 


2 Numerical Realization of a Photoacoustic In¬ 
version Formula 

Let U C be an open domain in and F a d — 1 dimensional hyperplane 
not intersecting U. Mathematically, photoacoustic imaging consists in solving 
the operator equation 

Q[/] 7*|rx{o,c3o) •> 

where / is a function with compact support in U and Q[/] is the trace on 
F X (0, oo) of the solution of the equation 




dttP — Ap = 0 in X (0, oo), 
p(-,0) = /(-) in 
dtp{-,Q) = 0 in . 

In other words, the photoacoustic imaging problem consists in identifying the 
initial source / from measurement data g = _p|rx(o,oo)- 

An explicit inversion formula for Q in terms of the Fourier transforms of / 
and g Q[/] has been found in [5S]. Let {x,y) S x K+. Assume without 
loss of generality (by choice of proper basis) that F is the hyperplane described 
by ?/ = 0. Then the reconstruction reads as follows: 

‘2K 

F[f] (K) = ^F[Q/] {K^, « (K)). (1) 

k(K) 

where F denotes the d-dimensional Fourier transform: 

r|/im 

and 

K{K)=s\gn{Ky)^Kl + Kl, 

K={K^,Ky) . 

Here, the variables x, are in whereas y, Ky G K. 

For the numerical realization these three steps have to be realized in discrete 
form: We denote evaluations of a function tp at sampling points {xm,yn) G 
(-X/2,X/2)^-i X (0,y) by 

^m,n •— Vn) ■ (2) 

For convenience, we will modify this notation in case of evaluations on an equi- 
spaced Cartesian grid. We dehne the d-dimensional grid 

G, X Gy := {-A,/2,..., N,/2 - x {0,..., - 1} , 

and assume our sampling points to be located on mA^juAy, where 

(m, 77.) G G 3 ; X Gy , 

and write 

— pipXlAxi TlAy) , (3) 

where A^, := X/N^ resp. Ay := Y/Ny are the occurring step sizes. 

In frequency domain, we have to sample symmetrically with respect to Ky. 
Therefore, we also introduce the interval 


Gk, := {-Ay/2,...,Ary/2-1}. 





Since we will have to deal with evaluations that are partially in-grid, partially 
not necessarily in-grid, we will also use combinations of and ([^. In this 
paper, we will make use of discretizations of the source function /, the data 
function g and their Fourier transforms / resp. g. 

Let in the following 


/bi = H 

{m,n)eGa:XGy 


^-2Tvi{j-m+ln)/{N^ Ny 


denote the d-dimensional discrete Fourier transform with respect to space and 
time. By discretizing formula Q via Riemann sums it follows 
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^—2tz\ Kj^n/Ny 


TlGGy 


27riij-m+ln)/Nt\ 

/ ^ ^ ym,n • 


Tn G Gx 


(4) 


where 


i^j,i = sign {l)\/p + P , 

{j,l) e G, X GKy ■ 

This is the formula from 

Remark 1. Note that we use the interval notation for the integer multi-indices 
for notational convenience. Moreover, we also choose the length of the Fourier 
transforms to be equal to in the first d — 1 dimensions, respectively. This 
could be generalised without changes in practice. 

Now, we assume to sample g at M , not necessarily uniform, points G 
(-X/2,X/2)‘^-b Then, 


nGG, 


— 27riKj jn/Ny 


M 




(5) 


The term hm represents the area of the detector surface around and has to 

M 

fulfil ^ hm = {Nx^xY~^ = ■ Note that the original formula Q can be 

ra—1 

received from by choosing {xm} to contain all points on the grid N^Gx. 

Formula (pT can be interpreted as follows: Once we have computed the 
Fourier transform of the data and evaluated the Fourier transform at non- 
equidistant points with respect to the third coordinate, we obtain the (stan¬ 
dard, equispaced) Fourier coefficients of /. The image can then be obtained by 
applying standard FFT techniques. 




The straightforward evaluation of the sums on the right hand side of ([^ 
would lead to a computational complexity of order Ny x Usually this is 
improved by the use of FFT methods, which have the drawback that they need 
both the data and evaluation grid to be equispaced in each coordinate. This 
means that if we want to compute (§ efficiently, we have to interpolate both in 
domain- and frequency space. A simple way of doing that is by using polynomial 
interpolation. It is used for photoacoustic reconstruction purposes for instance 
in the k-wave toolbox for Matlab m- Unfortunately, this kind of interpolation 
seems to be sub-optimal for Fourier-interpolation with respect to both accuracy 
and computational costs [3 US] 

A regularized inverse k-space interpolation has already been shown to yield 
better reconstruction results m- The superiority of applying the NUFFT, com¬ 
pared to linear interpolation, has been shown theoretically and computationally 

by i. 


3 The non-uniform fast Fourier transform 


This section is devoted to the brief explanation of the theory and the applica¬ 
bility of the non-uniform Fourier transform, where we explain both the NER- 
NUFFT (subsection 3.1) and the NED-NUFFT (subsection |3.2| ) in the form 
(and spatial dimensions) we utilise them afterwards. 

The NEDNER-NUFFT algorithm used for implementing (§ essentially (up 
to scaling factors) consists of the following steps: 


1. Compute a d — 1 dimensional NED-NUFFT in the a;-coordinates due to 
our detector placement. 

2. Compute a one-dimensional NER-NUFFT in the ATy-coordinate as indi¬ 
cated by the reconstruction formula (§. 

3. Compute an equispaced d-dim inverse FFT to obtain a d dimensional 
picture of the initial pressure distribution. 


3.1 The non-equispaced range (NER-NUFFT) case 

With the NER-NUFFT (non equispaced range - non-uniform FFT) it is possible 
to efficiently evaluate the discrete Fourier transform at non-equispaced positions 
in frequency domain. 

To this end, we introduce the one dimensional discrete Fourier transform, 
evaluated at non-equispaced grid points ki G M: 

Z = 1,...,M. (6) 

n^Gy 

In order to find an efficient algorithm for evaluation of ([^, we use a window 
function dt, an oversampling factor c > 1 and a parameter c < a < tt{2c — 1) 
that satisfy: 



1. 'I' is continuous inside some finite interval [—a, a] and has its support in 
this interval and 


2. 'I' is positive in the interval [—7r,7r]. 

Then (see di) we have the following representation for the Fourier modes 
occurring in ([^: 




c 




\9\ < TT . 


(7) 


By assumption, both and are concentrated around 0. So we approximate 
the sum over all A: G Z by the sum over the 2K integers k that are closest to 
Ki + k. By choosing 9 = 2'Kn/N — tt and inserting Q in ([^, we obtain 


« E E 

k=-K+l nGGy " 

I = . 


( 8 ) 


Here K denotes the interpolation length and 
:= ^'(27rn/iVy - tt) , 

V 27r 

where is the nearest integer (i.e. the nearest equispaced grid point) to Ki + k. 

The choice of is made in accordance with the assumptions above, so we 
need 'I' to have compact support. Furthermore, to make the approximation 
in (j^ reasonable, its Fourier transform d' needs to be concentrated as much 
as possible in [—K,K]. In practice, a common choice for is the Kaiser- 
Bessel function, which fulfils the needed conditions, and its Fourier transform is 
analytically computable. 


3.2 The non-equispaced data (NED-NUFFT) case 

A second major aim of the present work is to handle data measured at non- 
equispaced acquisition points in an efficient and accurate way. Therefore we 
introduce the non-equispaced data, d — 1 dimensional DFT 

M 

j G Gx • 


The theory for the NED-NUFFT is largely analogous to the NER-NUFFT [7] as 
described in Subsection 3.1 The representation 0 is here used for each entry 





of j and inserted (with now setting 9 = 2?™ /TV) into formula , which leads 
to 



■' m=l 


( 10 ) 


where the entries in firn,k are the nearest integers to Xm + k. Here we have used 
the abbreviations 


d-l 


:= n , 



for the needed evaluations of and 

Further remarks on the implementation of the NED- and NER-NUFFT, as 
well as a summery about the properties of the Kaiser-Bessel function and its 
Fourier transform can be found in [TUg. 


4 Comparison of NER-NUFFT and k-wave FFT 


Before we turn to the evaluation of the algorithm we describe the photoacoustic 
setup. Our device consists of a FP (Fabry Perot) polymer film sensor for inter¬ 
rogation mm- A 50 Hz pulsed laser source and a subsequent optical parametric 
oscillator (OPO) provide optical pulses. These pulses have a very narrow band¬ 
width and can be tuned within the visible and near infrared spectrum. The 
optical pulses are then transmitted via an optical fibre. When the light is emit¬ 
ted it diverges and impinges upon a sample with homogeneous fluence, thus 
generating a photoacoustic signal. This signal is then recorded via the FP- 
sensor head. The sensor head consists of an approximately 38 pm thick poly¬ 
mer (Parylene C) which is sandwiched between two dichroic dielectric coatings. 
These dichroic mirrors have a noteworthy transmission characteristic. Light 
from 600 to 1200 nm can pass the mirrors largely unabated, whereas the reflec¬ 
tivity from 1500 to 1650 nm (sensor interrogation band) is about 95% [37]. The 
incident photoacoustic wave produces a linear change in the optical thickness of 
the polymer film. A focused continuous wave laser, operating within the inter¬ 
rogation band, can now determine the change of thickness at the interrogation 
point via FP-interferometry. 

We choose two 2D targets for comparison, a star and a USAF (US Air Force) 
resolution test chart. Both targets are made of glass with a vacuum-deposited 
durable chromium coating. The star target has 72 sectors on a pattern diameter 
of 5 mm, with an unresolved core diameter of 100 pm. 
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Figure 1: Segment of the MIP (maximum intensity projection) in the z-axis of 
a star sample, reconstructed using FFT and linear interpolation (top left) and 
using the NUFFT (top right) with a square plotted around the center. The 
intensity of the reconstructed image along the sidelines of the square is plotted, 
for both reconstructions. The purple line indicates the frequency of the star 
sample. 


The targets are positioned in parallel to the sensor surface at a distance of 
about 4 mm, and water is used as coupling medium between the target and the 
sensor. 

The interpolation length for the NER-NUFFT reconstruction is K = 6. The 
computational times are shown in Tab. [^showing that the linearly interpolated 
FFT is about 30 % slower than the NER-NUFFT. 

A segment of the reconstructed star target is shown in Fig. The intensity 
is plotted against the sides of an imaginary square (2.67 mm) placed around the 
center of the star phantom. It is clearly visible, that the FFT reconstruction 
is not able to represent the line pairs, when the density exceeds lOlp/mm, 
corresponding to a resolution of 100 pm, whereas they are still largely visible for 
the NER-NUFFT reconstruction. 

For a quantitative comparison of the resolution we use the USAF chart. It is 
recorded on an area of 146 x 146 sensor points corresponding to 1.022 x 1.022 cm^ 
with a grid spacing of 70 pm and a time resolution of 8 ns. As yet there is no 
standardized procedure to measure the resolution of a photoacoustic imaging 



















Table 1: Comparison between the NED-NUFFT and FFT reconstruction for a 
USAF-chart and a comparison of computational times for both phantoms. The 
improvement in percent was caluclated by: 100 x (1 — FFT/NUFFT) 



NUFFT 

FFT 

Improvement 

FWHM axial LSF 

23.23 ± 0.56 pm 

62.34 ± 0.62 pm 

168.47 ±6.88% 

FWHM lateral LSF 

33.44 ± 7.95 pm 

40.82 ± 7.34 pm 

18.63 ±8.50% 

Time: Star target 

140 s 

189 s 

35% 

Time: USAF chart 

298 s 

384 s 

29% 


system. We proceed similar to m, by fitting a line spread function (LSF) and 
an edge spread function (FSF) to the intensities of our reconstructed data. For 
the LSF to be meaningful, its source has to approximate a spatial delta function. 
This is the case in the z-axis, since the chrome coating of the USAF target is 
just about 0.1 p.m thick. 

We fit a Cauchy-Lorentz distribution to the z-axis of our reconstructed data 
for 49 adjacent a; 2 /-coordinates. Their positions are marked as a black square 
within the white square, depicted in the bottom images of Fig. The recon¬ 
struction in the z-axis for a single point is shown in the bottom right inlay of 
Fig. [21 for 8 points, covering a distance of 94.74 pm in the z-direction, around 
the maximum intensity. A fit of the Lorentz distribution is shown for both re¬ 
construction methods. The FWHM {full width, have maximum) of the Lorentz 
distribution. 


I{z) 


2aow 

TT -I- 4(z — Zq)^) ’ 


is the parameter w. The output /(z) is the intensity in dependence from the 
z-axis, and Zq and Oq are fitting parameters. 

The average and standard deviation of w for both of the 49 datasets are 
shown in Table The line spread function FWHM of the FFT-reconstruction 
turns out to be more than twice as big as the one of the NFR-NUFFT recon¬ 
struction. 

For the lateral resolution, there is no target approximating a delta function, 
so we have to use the FSF instead: 


I{x) = lo+ao 


fl 

— arctan 
Vtt 


X — Xq 
wl2 



The w here is the FWHM of the associated LSF, and Jq, a:o and oq are fitting 
parameters. The ESF requires a step function as source, of which our target 
provides plenty. We choose the long side of 12 bars, marked by white lines 
in the bottom images of Fig. for this fit. The data for a particular line 
is shown in the bottom left inlay. We omitted all datasets, where only one 
point marked the transition from low to high intensity, rendering the edge fit 
unreliable and resulting in unrealistic improvements of our new method well 

















Figure 2: Segment of the xy-MIP of an USAF chart reconstruction conducted 
with FFT (left) and NER-NUFFT (right). On the bottom images the data-sets 
used for the quantitative resolution are marked. The black square depicts the 
49 xy-coordinates used for the LSF fit along the z-axis for the axial resolution. 
The data for a single point is shown in the bottom right inlay. The white lines 
show the intensity fit for the lateral resolution. The bottom left inlay shows the 
data for a single white line. 











over 100%. Finally we averaged over 15 edges. The results are shown in Table 
While the deviation between different edges is quite high, the improvement 
of our NER-NUFFT reconstruction for the 15 evaluated edges ranged from 6 to 
44%. 


5 Non-equidistant grid sampling 

The current setups allow data acquisition at just one single sensor point for 
each laser pulse excitation. Since our laser is operating at 50 Hz data recording 
of a typical sample requires several minutes. Reducing this acquisition time 
is a crucial step in advancing Photoacoustics towards clinical and preclinical 
application. Therefore, in this work we try to maximize the image quality for 
a given number of acquisition points. We are able to do this, because our 
newly implemented NEDNER-NUFFT is ideal for dealing with non-equispaced 
positioned sensors. This newly gained flexibility of sensor positioning offers 
many possibilities to enhance the image quality, compared against a rectangular 
grid. Eor instance a hexagonal grid was found to yield an efficiency of 90.8% 
compared with 78.5% for an exact reconstruction of a wave-number limited 
function m- 

Also any non-equispaced grids that may arrive from a specific experimental 
setup can be efficiently computed via the NEDNER-NUEET approach. 

Here, we will use it to tackle the limited view problem. Many papers deal 
with the limited view problem, when reconstructing images [6l|26l[l6]. Our 
approach to deal with this problem is different. It takes into account that in 
many cases the limiting factor is the number of sensor points and the limited 
view largely a consequence of this limitation. 

In our approach we therefore use a grid arrangement that is dense close to a 
center of interest and becomes sparser the further away the sampling points are 
located. We realize this by means of an equiangular, or equi-steradian sensor 
arrangement, where for a given point of interest each unit angle or steradian 
gets assigned one sensor point. 

To understand the limited view problem, it is helpful to define a detection 
region. According to [3^, this is the region which is enclosed by the normal 
lines from the edges of the sensor. Pressure waves always travel in the direction 
of the normal vector of the boundary of the expanding object. Therefore certain 
boundaries are invisible to the detector as depicted in Eig. 

5.1 Equiangular and equi-steradian projection sensor mask 

For the equiangular sensor arrangement a point of interest is chosen. Each line, 
connecting a sensor point with the point of interest, encloses a fixed angle to its 
adjacent line. In that sense we mimic a circular sensor array on a straight line. 
The position of the sensor points can be seen on top of the third image in Eig. 

01 



Figure 3: Depiction of the limited view problem. Edges whose normal vector 
cannot intersect with the sensor surface are invisible to the sensor. The invisible 
edges are the coarsely dotted lines. The detection region is marked by a grey 
background. The finely dotted lines are used to construct the invisible edges. 
Edges perpendicular to the sensor surface are always invisible for a plane sensor. 


The obvious expansion of an equiangular projection to 3D is an equi-steradian 
projection. This problem is analogous to the problem of placing equispaced 
points on a 3D sphere and then projecting the points, from the center of the 
sphere, through the points, onto a 2D plane outside the sphere. 

We developed an algorithm for this problem, which is explained in detail in 
Appendix Our input variables are the grid size, the distance of the center 
of interest from the sensor plane and the desired number of acquisition points, 
which will be rounded to the next convenient value. 

A sensor arrangement with 1625 points on a 226 x 226 grid is shown in the 
top left image in Eig. 

5.2 Weighting term 

To determine the weighting term in Eq. for 3D we introduce a function 
that describes the density of equidistant points per unit area Pp. In our specific 
case, Pp describes the density on a sphere around a center of interest. Eurther 
we assume that pp is spherically symmetric and decreases quadratically with 
the distance from the center of interest r: pp^s oc 1/r^. We now define Pp^m for 
a plane positioned at distance tq from the center of interest. In this case Pp^sif) 
attenuates by a factor of sin a, where a = arcsin(ro/r) is the angle of incidence. 
Hence Pp^m cx By applying the spacing of the regular grid . This yields 

a weighting term of: 

hm{r) oc 

Analogously we can derive for 2D: 

hm{r) oc 

We applied a normalization after the reconstruction to all measurements. 
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Figure 4: Various reconstructions of a tree phantom (top) with different sensor 
arrangements. All sensor arrangements are confined to 32 sensor points. The 
sensor positions are indicated as white rectangles on the top of the images. 
The second image shows the best (see Figj^ equispaced sensor arrangement, 
with a distance of 13 points between each sensor. The third image shows the 
NEDNER-NUFFT reconstruction with equiangular arranged sensor positions. 
The bottom image shows the same sensor arrangement, but all omitted sensor 
points are linearly interpolated and afterwards a NER-NUFFT reconstruction 
was conducted. 


For the application of this method to the FP setup it is noteworthy that 
there is a frequency dependency on sensitivity which itself depends on the angle 
of incidence, which has been extensively discussed in [^. The angle of inci¬ 
dence for our specific setup is 62°. At this angle, the frequency components 
around 2 MHz get attenuated by more than 10 dB. Below IMHz the frequency 
response remains quite stable (attenuation below 5 dB) for the measurement 
angles occurring in our setup. 






6 Computational assessment of different non- 
equispaced grid arrangements in two dimen¬ 
sions to tackle the limited view problem 

A tree phantom, designed by Brian Hurshman and licensed under CC BY 3.(Q is 
chosen for the 2 dimensional computational experiments on a grid with x = 1024 
z = 256 points. A forward simulation is conducted via k-wave m- The forward 
simulation of the k-wave toolbox is based on a first order k-space model. A PML 
(perfectly matched layer) of 64 gridpoints is added, as well as 30 dB of noise. 

In Fig. our computational phantom is shown at the top. For each recon¬ 
struction a subset of 32 out of the 1024 possible sensor positions was chosen. In 
Fig. i their positions are marked at the top of each reconstructed image. For 
the equispaced sensor arrangements, we let the distance between two adjacent 
sensor points sweep from 1 to 32. The sensor points where always centered in 
the x-axis. 

To compare the different reconstruction methods we used the correlation 
coefficient and the Tenenbaum sharpness. These quality measures are explained 
in Appendix |10| 

We applied the correlation coefficient only within the region of interest 
marked by the white circle in Fig. The Tenenbaum sharpness was calcu¬ 
lated on the smallest rectangle, containing all pixels within the circle. The 
results are shown in Fig. 

The Tenenbaum sharpness for the equiangular sensor placement was 23001, 
which is above all values for the equispaced arrangements. The correlation 
coefficient was 0.913 compared to 0.849, for the best equispaced arrangement. 
In other words, the equiangular arrangement is 42.3 % closer to a full correlation, 
than any equispaced grid. 

In Fig. I^the competing reconstructions are compared. While the crown of 
the tree is depicted quite well for the equispaced reconstruction, the trunk of the 
tree is barely visible. This is due to the limited view problem. When the equis¬ 
paced interval increases, the tree becomes visible, but at the cost of the crown’s 
quality. In the equiangular arrangement a trade off between these two effects is 
achieved. Additionally the weighting term for the outmost sensors is 17 times 
the weighting term for the sensor point closest to the middle. This amplifies the 
occurrence of artefacts, particularly outside of our region of interest. 

The bottom image in Fig. shows the equiangular sensor arrangement, 
reconstructed in a conventional manner. The missing sensor points are inter¬ 
polated to an equispaced grid, and a NER-NUFFT reconstruction is applied 
afterwards. We conducted a linear interpolation from our subset to all 1024 
sensor points. The correlation coefficient for this outcome was 0.7348 while 
the sharpness measure was 21474. This outcome exemplifies the clear supe¬ 
riority of the NUFFT to conventional FFT reconstruction when dealing with 
non-equispaced grids. 


^ http: //thenounproject. com/term / tree /16622/ 
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Figure 5: Correlation coefRcient and Tenenbaum sharpness for equispaced sen¬ 
sor arrangements with intervals between the sensor points reaching from 1 to 
32. The maximum of the correlation coefficient is at 13. The corresponding 
reconstruction is shown in Fig. The straight lines indicate the results for the 
equiangular projection. 


7 3D application of the NED-NER-NUFFT with 
real data 

For a qualitative assessment of our new sensor arrangement we need a 3D phan¬ 
tom. We choose a yarn which we record on a rectangular grid with 226 x 226 = 
51076 sensor points, with a grid spacing of 60 pm and time sampling of dt = 8 ns. 
Hence an area of 13.56 x 13.56 mm is covered. As coupling medium water is 
used, in which the yarn is fully immersed. 

To determine the utility of non-equispaced grid sampling, we follow a certain 
routine. First we acquire a densely sampled dataset. Then we use a very small 
subset of the initially collected sensor data, to test different sensor arrangements. 
Therefore we can always use the complete reconstruction as our model standard 
together with the quality measurements explained in Appendix |10| 

An upsampling factor of 2 was used for all reconstructions, hence the recon¬ 
structed image of the MIP for the xy plane consists of 452 x 452 pixels. The 
complete reconstruction with the NER-NUFFT took 154 seconds. Maximum 
intensity projections (MIP) of this full reconstruction for all axis are shown in 
Fig. § The MIP in the xy plane is our model standard, for comparison with 
the other reconstructions. 

For the equi-steradian sensor mask we choose our center of interest right 
in the center of the xy-MIP where the little knot can be seen, 3.6 mm off the 
sensor surface. The resulting sensor mask, including the reconstruction is shown 
in Fig. it consists of 1625 sensor points (or 3.18 % of the initial number of 
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Figure 6: MIPs of all 3 planes for the NER-NUFFT reconstruction of a yarn 
phantom. 


sensor points). The weighting term, accounting for sensor sparsity is 9.7 times 
higher for the outermost sensor point, than for centermost sensor points in the 
a:?/-plane. The reconstruction for this arrangement took 134 seconds. 

We compared this arrangement to rectangular grids, which all had 41 x 41 = 
1681 sensor points (or 3.29 % of the initial number of sensor points) but varying 
distances between two adjacent points. The grid with an interval between sensor 
points of 5 is shown on the top left in Fig. 

In Fig. [^the equi-steradian grid is compared with 2 equispaced grids. To 
get a more precise measure of the correlation between the reconstructed images 
and the model standard, we calculated the correlation coefficient only within a 
region of interest. The region of interest is firstly confined by a centered disc, 
whose boundary is shown as a dotted circle in the inlay in Fig. The a:-axis 
shows the number of pixels of interest used to calculate the correlation coefficient 
within this disc. These pixels are increased according to the intensity of the 
corresponding pixels of the model’s standard MIP. Fig. [^demonstrates that the 
correlation coefficient for the equi-steradian arrangement, always remains better, 
within the depicted disc, when competing against the two strongest equispaced 
grids. 

In Fig. I^the 4 equispaced arrangements, with intervals reaching from 2 to 5 
are compared to the equi-steradian grid. Here the pixels of interest are set to a 
threshold of at least 4% of the maximum value, while the diameter is increased. 
The dotted straight lines indicate the side length of the square of a particular 



Figure 7: The sensor placement for the cirular arrangement is shown on the 
top left image, comprising 1625 sensor points. On the top right an equispaced 
sensor arrangement with 41 x 41 = 1689 sensor points is displayed. The intervall 
between two adjacent sensor points is 5 for this configuration. The bottom im¬ 
ages show MIPs of the NED-NER-NUFFT reconstruction only using the sensor 
points shown above. 



pixels used for correlation coefficient calculation 


Figure 8: Correlation coefficient calculated on different number of pixels, for 
3 different reconstructed MIPs. Pixels of interest are added according to the 
intensity of the corresponding pixels of the model standard image (Fig. [^. 
In the inlay the correlation coefficient mask is shown for 20087 points which 
corresponds to all pixels with at least 4% of the maximum value of the model 
standard and is the value used in Fig. 1^ 



















Figure 9: Correlation Coefficient calculated on pixels of interest with at least 
4% of the maximum value of the model standard, for an increasing diameter of 
the confining disc. The straight line at 8.4 mm corresponds to the straight line 
in Fig. 1^ The other straight lines indicate the side length of the square of a 
particular equispaced sensor point arrangement. 


equispaced sensor point arrangement. As expected the correlation coefficient 
for this equispaced sensor arrangement starts to decline around that threshold. 

The correlation coefficient for the equi-steradian grid starts to fall behind 
the equispaced grid of interval 5 towards the end. This is expected, since the 
equi-steradian grid is only meant to give better results for a region of interest 
around the center. This is very clearly the case. Between a diameter of 2 to 
8 mm the correlation coefficient is on average 25.14 % closer to a value of 1 
than its strongest equispaced contender of interval 4, and 30.87 % better than 
the interval 5 grid. At a diameter of 7.2 mm the correlation coefficient for the 
equi-steradian grid was 0.960 compared to 0.947 for the interval 4 grid and 0.944 
for the interval 5 grid. This is 24.8 % and 28.6% closer to full correlation. 

Fig. [^shows the normalized Tenenbaum sharpness. The Tenenbaum sharp¬ 
ness, unlike the correlation coefficient, cannot be calculated on non-adjacent grid 
points, therefore it has been calculated on the smallest rectangle, that contains 
all pixels of interest. The equi-steradian has the highest sharpness for most 
values, with the interval 4 grid being very slightly better around a diameter of 
3-4 mm. There is a drop of the Tenenbaum sharpness towards the end. 


8 Discussion and Conclusion 

We computationally implemented a 3D non-uniform FFT photoacoustic image 
reconstruction, called NER-NUFFT (non equispaced range-non uniform FFT) 
to efficiently deal with the non-equispaced Fourier transform evaluations arising 
in the reconstruction formula. This method was compared with the k-wave 
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Figure 10: Tenenbaum Sharpness calculated on pixels of interest with at least 
4% of the maximum value of the model standard. The Tenenbaum sharpness is 
calculated on the smallest rectangle, that contains all pixels of interest. 


implemented FFT reconstruction, which uses a polynomial interpolation. The 
two reconstruction methods where compared using 2D targets. The lateral 
resolution showed an improvement of 18.63 ± 8.5% which is in good agreement 
with the illustrative results for the star target. The axial resolution showed an 
improvement of 168.47 ± 6.88%. The computation time was about 30% less, 
for the NER-NUFFT than the linearly interpolated FFT reconstruction. In 
conclusion the NER-NUFFT reconstruction proved to be unequivocally superior 
to conventional linear interpolation FFT reconstruction methods. 

We further implemented the NED-NER-NUFFT (non equispaced data-NER- 
NUFFT), which allowed us to efficiently reconstruct from data recorded at non- 
equispaced placed sensor points. This newly gained flexibility was used to tackle 
the limited view problem, by placing sensors more sparsely further away from 
the center of interest. We developed an equiangular sensor placement for 2D 
and an equi-steradian placement in 3D, which assigns one sensor point to each 
angle/steradian for a given center of interest. In the 2D computational simula¬ 
tions we showed that this arrangement significantly enhances image quality in 
comparison to regular grids. 

In 3D we conducted experiments, where a yarn phantom was recorded. The 
maximum intensity projection (MIP) of the full reconstruction was compared to 
MIPs of reconstructions that only used about 3% of the original data. Within 
our region of interest, the correlation of our image was 0.96, which is 24.8% 
closer to full correlation than the best equispaced arrangement, reconstructed 
from slightly more sensor points. 

The sensor placement to tackle the limited view problem, combined with 
the NED-NER-NUFFT gives significantly better results for an object located 
at the center of a bigger sensor surface. This result was confirmed in the 2D 
simulation as well as for real data in 3D. 












Acknowledgment 

This work is supported by the Medical University of Vienna, the European 
projects FAMOS (FP7 ICT 317744) and FUN OCT (FP7 HEALTH 201880), 
Macular Vision Research Foundation (MVRF, USA), Austrian Science Fund 
(FWF), Project P26687-N25 (Interdisciplinary Coupled Physics Imaging), and 
the Christian Doppler Society (Christian Doppler Laboratory ’’Laser develop¬ 
ment and their application in medicine”). 


References 

[1] P. Beard. Biomedical photoacoustic imaging. Interface Focus, 1:602-631, 
2011 . 

[2] P.C. Beard. Two-dimensional ultrasound receive array using an angle-tuned 
fabry-perot polymer Him sensor for transducer field characterization and 
transmission ultrasound imaging. IEEE Trans. Ultrason., Eerroeletr., Freq. 
Control, 52(6):1002-1012, June 2005. 

[3] P.C. Beard, F. Perennes, and T.N. Mills. Transduction mechanisms of the 
fabry-perot polymer film sensing concept for wideband ultrasound detec¬ 
tion. IEEE Trans. Ultrason., Ferroeletr., Freq. Control, 46(6):1575-1582, 
November 1999. 

[4] P. Burgholzer, C. Hofer, G. Paltauf, M. Haltmeier, and O. Scherzer. 
Thermoacoustic tomography with integrating area and line detectors. IEEE 
Trans. Ultrason., Eerroeletr., Ereq. Control, 52(9):1577-1583, September 
2005. 

[5] B.T. Cox and P.C. Beard. The frequency-dependent directivity of a pla¬ 
nar fabry-perot polymer film ultrasound sensor. IEEE Trans. Ultrason., 
Ferroeletr., Freq. Control, 54(2):394-404, 2007. 

[6] W. Dan, C. Tao, X.-J. Liu, and X.-D. Wang. Influence of limited-view 
scanning on depth imaging of photoacoustic tomography. Chinese Physics 
B, 21(1):014301, 2012. 

[7] K. Fourmont. Non-equispaced fast Fourier transforms with applications to 
tomography. J. Fourier Anal. AppL, 9(5):431-450, 2003. 

[8] F. C. A. Groen, 1. T. Young, and G. Ligthart. A comparison of different 
focus functions for use in autofocus algorithms. Cytometry, 6(2):81-91, 
1985. 

[9] M. Haltmeier, O. Scherzer, and G. Zangerl. A reconstruction algorithm for 
photoacoustic imaging based on the nonuniform FFT. IEEE Trans. Med. 
Imag., 28(11):1727-1735, November 2009. 



[10] M. Jaeger, S. Schiipbach, A. Gertsch, M. Kitz, and M. Frenz. Fourier 
reconstruction in optoacoustic imaging using truncated regularized inverse 
k-space interpolation. Inverse ProbL, 23:S51-S63, 2007. 

[11] P. Kuchment and L. Kunyansky. Mathematics of thermoacoustic tomogra¬ 
phy. European J. Appl. Math., 19:191-224, 2008. 

[12] G. Paltauf, R. Nuster, M. Haltmeier, and P. Burgholzer. Photoacoustic 
tomography with integrating area and line detectors. In L. V. Wang, editor. 
Photoacoustic Imaging and Spectroscopy, Optical Science and Engineering, 
pages 251-263. GRC Press, Boca Raton, FL, 2009. 

[13] D. P. Petersen and D. Middleton. Sampling and reconstruction of wave¬ 
number-limited functions in Wdimensional Euclidean spaces. Inform, and 
Control, 5:279-323, 1962. 

[14] O. Scherzer, editor. Handbook of Mathematical Methods in Imaging. 
Springer, New York, 2011. 

[15] R. Schulze, G. Zangerl, M. Holotta, D. Meyer, F. Handle, R. Nuster, G. Pal¬ 
tauf, and O. Scherzer. On the use of frequency-domain reconstruction al¬ 
gorithms for photoacoustic imaging. J. Biomed. Opt., 16(8):086002, 2011. 
Funded by the Austrian Science Fund (FWF) within the FSP S105 - “Pho¬ 
toacoustic Imaging”. 

[16] G. Tao and X. Liu. Reconstruction of high quality photoacoustic tomog¬ 
raphy with a limited-view scanning. Opt. Express, 18(3):2760-2766, Feb 
2010 . 

[17] B. E. Treeby and B. T. Gox. K-Wave: MATLAB toolbox for the simu¬ 
lation and reconstruction of photoacoustic wace fields. J. Biomed. Opt., 
15:021314, 2010. 

[18] B. E. Treeby, T. K. Varslot, E. Z. Zhang, J. G. Laufer, and P. C. Beard. 
Automatic sound speed selection in photoacoustic image reconstruction us¬ 
ing an autofocus approach. Biomed. Opt. Express, 16(9):090501-090501-3, 
2011 . 

[19] L. V. Wang, editor. Photoacoustic Imaging and Spectroscopy. Optical 
Science and Engineering. GRG Press, Boca Raton, 2009. 

[20] M. Xu and L. V. Wang. Exact frequency-domain reconstruction for ther¬ 
moacoustic tomography-I: Planar geometry. IEEE Trans. Med. Imag., 
21:823-828, 2002. 

[21] M. Xu and L. V. Wang. Time-domain reconstruction for thermoacoustic 
tomography in a spherical geometry. IEEE Trans. Med. Imag., 21(7):814- 
822, 2002. 



[22] M. Xu and L. V. Wang. Analytic explanation of spatial resolution related to 
bandwidth and detector aperture size in thermoacoustic or photoacoustic 
reconstruction. Phys. Rev. E, 67(5):0566051-05660515, 2003. 

[23] M. Xu and L. V. Wang. Universal back-projection algorithm for photoa¬ 
coustic computed tomography. Phys. Rev. E, 71(1), 2005. 

[24] M. Xu and L. V. Wang. Photoacoustic imaging in biomedicine. Rev. Sci. 
Instruments, 77(4), 2006. 

[25] Y. Xu, D. Feng, and L. V. Wang. Exact frequency-domain reconstrcution 
for thermoacoustic tomography — 1: Planar geometry. IEEE Trans. Med. 
Imag., 21(7):823-828, 2002. 

[26] Y. Xu, L. V. Wang, G. Ambartsoumian, and P. Kuchment. Reconstructions 
in limited-view thermoacoustic tomography. Med. Phys., 31(4):724-733, 
2004. 

[27] E. Zhang, J. Laufer, and P. Beard. Backward-mode multiwavelength pho- 
toacoustic scanner using a planar Eabry-Perot polymer film ultrasound 
sensor for high-resolution three-dimensional imaging of biological tissues. 
App. Opt, 47:561-577, 2008. 

9 Appendix: Algorithm for equi-steradian sen¬ 
sor arrangement 

In our algorithm, the grid size and the distance of the center of interest from 
the sensor plane is defined. The number of sensor points N will be rounded to 
the next convenient value. 

Our point of interest is placed at z = r^, centered at a square xy grid. The 
point of interest is the center of a spherical coordinate system, with the polar 
angle 0 = 0 at the 2 :-axis towards the xy-grid. 

Eirst we determine the steradian 12 of the spherical cap from the point of 
interest, that projects onto the acquistion point plane via 

12 = 27r(l - cos (Omax)) ■ 

This leads to a unit steradian oj = tl/N with N being the number of sensors 
one would like to record the signal with. The sphere cap is then subdivided into 
slices k which satisfy the condition 

w jfc = 27r(cos(6»fc-i) - cos(6»fc)) , 

where 9i encloses exactly one unit steradian oj and jk has to be a power of two, in 
order to guarantee some symmetry. The value of jk doubles, when > 1.8 • r^, 
where is the chord length between two points on k and is the distance to 
the closest point on fc — 1. These values are chosen in order to guarantee We 



apply some restrictions to approximate equidistance between acquisition points 
on the sensor surface. 

The azimuthal angles for a slice k are calulated according to: 

= (27ri) /jfc + 7r/jfc + , 

with i = 0,..., j — 1, where 

+ (fc - 1) 27r/(jfc_i) 

stems from the former slice k — 1 . The sensor points are now placed on the 
a:?/-plane at the position indicated by the spherical angular coordinates: 


(pol, az) = {{Ok + Ok+i) /2, pi^k) 


10 Appendix: Quality measures 

The correlation coefficient p is a measure of the linear dependence between two 
images Ui and U 2 - 

Its range is [—1,1] and a correlation coefficient close to 1 indicates linear 
dependence [T3]. It is defined via the variance, Var([/j) of each image and the 
covariance, Cov (I/i, U 2 ) of the two images: 



( 11 ) 


We did not choose the widely used distance measure because it is a 
morphological distance measure, meaning it defines the distance between two 
images by the distance between their level sets. Therefore two identical linearly 
dependent images can have a correlation coefficient of 1 and still a huge 
distance. Normalizing the images only mitigates this problem, because single, 
high intensity artifacts or a varying intensity over the whole image can greatly 
alter the distance. 

While the correlation coefficient is a good measure for the overall similarity 
between two images it does not include any sharpness measure. Hence blurred 
edges are punished very little, in comparison to slight variations of homogeneous 
areas. To address this shortcoming we chose a sharpness measure or focus 
function as a second quality criterion. Sharpness measures are obtained from 
some measure of the high frequency content of an image [5]. They have also 
been used to select the best sound speed in photoacoustic image reconstruction 
[H] . Out of the plethora of published focus functions we select the Tenenbaum 
function, because of its robustness to noise: 



( 12 ) 


with g as the Sobel operator: 



( 13 ) 




Like the norm and unlike the correlation coefficient the Tenenbaum 
function is an extensive measure, meaning it increases with image dimensions. 
Therefore we normalized it to J^Tenenbaum = ^Tenenbaum/.^, where N is the 
number of elements in U. 



